%--------------------------------------------------------------------------
% function plots the neutral stability curves and then, for a given set of
% parameter values (k, Ma), it traces out the path in 'stability space'
% that results from the time evolution of the problem.
%
% can also compute the eigenvalues along this path
%--------------------------------------------------------------------------

function make_path(k0, N_eigs)

dispersion = 1;

close all

if nargin == 1
N_eigs = 0;
end

addpath ../
p = params;

% k = linspace(0, 3, 100);
k = logspace(-2, log10(3), 100);

% compute the dispersion relation
if dispersion
    subplot(1,2,1);
    ev = comp_eigs(k, p, 1);
    
    plot(k, ev, 'k', 'linewidth',1);
    hold on;
    plot([k(1), k(end)], [0 0], 'k');
    
    
%     semilogx(k, ev, 'k', 'linewidth',1);
%     hold on;
%     semilogx([k(1), k(end)], [0 0], 'k');
    
    
    
    xlabel('$\check{k}$','interpreter','latex','fontsize',12);
    ylabel('$\check{\lambda}$','interpreter','latex','fontsize',12);
    
    subplot(1,2,2);
end

% plot the neutral stability curves
[k, Ma] = exact_neutral_stab(k, 0, 0);
plot(k, Ma * p.delta * p.beta * (1 - p.beta), 'k','linewidth',1);


ylim([0, 130]);
xlabel('$\check{k}$','interpreter','latex','fontsize',12);
ylabel('$\mathcal{M}$', 'interpreter','latex','fontsize',12);

if dispersion
    hold on;
    plot([k(1) k(end)], [p.Ma_prime, p.Ma_prime],'k-.');
    
%     semilogx([k(1) k(end)], [p.Ma_prime, p.Ma_prime],'k-.');
end
   

% compute h
t = [logspace(-6, -2, 20) linspace(1.1e-2, p.Tmax, 100)];
% t = logspace(-6, log10(p.Tmax), 100);
h = -(1 + lambertw(-p.beta / (p.beta - 1) * exp(((-p.beta + p.delta * t) / (p.beta - 1))))) * (p.beta - 1);

% compute time-dependent wavenumber and Marangoni number

Ma = (h - 1 + p.beta) / p.beta * p.Ma;
k = k0 * h;

hold on;
plot(k, p.delta * p.beta * (1 - p.beta) * Ma, 'k--','linewidth',1);
plot(k(1), p.delta * p.beta * (1 - p.beta) * Ma(1), 'k*');


% compute the eigenvalues on this path
if N_eigs > 0
    ev = zeros(N_eigs, length(t));
    
    for i = 1:length(t)
        
%         p.Ma = Ma(i);
        [tmp1, tmp2, L] = comp_eigs(k0, p, h(i));
        if (i == 1)
            [ef, tmp1] = eigs(L, 1, 'SM');
        end
           
        ev(:,i) = eigs(L, N_eigs, 'SM');
    end
    
    % plot the eigenvalues
    figure;
    plot(t, ev);
    xlabel('$t$','interpreter','latex','fontsize',12);
    ylabel('$\check{\lambda}$','interpreter','latex','fontsize',12);
    
    
    % compute the amplificiation factors
    lambda = ev(1,:);
    
    li = @(X) interp1(t, lambda, X);
    A_asy_i = @(T) exp(quad( @(t) li(t), t(1), T));
    
        
    lin = stab(k0, p, ef / norm(ef, 'inf'));
    A = max(abs(lin.y), [], 1);
    
    sig = diff(log(A)) ./ diff(lin.x);
    
    
    % plot the top eigenvalue
    
    ti = linspace(t(1), p.Tmax, 30);
    til = logspace(log10(t(1)), log10(p.Tmax), 30);
    figure
    subplot(2,1,1);
    semilogx(lin.x(1:end-1), sig, 'k'); hold on;
    semilogx(til, interp1(t, lambda, til), 'k*')
    xlabel('$t$');
    ylabel('Top eigenvalue');
    subplot(2,1,2);
    plot(lin.x(1:end-1), sig, 'k'); hold on;
    plot(ti, interp1(t, lambda, ti), 'k*')
    xlabel('$t$');
    ylabel('Top eigenvalue');
    
    figure;
    plot(lin.x, A, 'k','linewidth',1);
    hold on;
    
    for i = 1:length(ti)
        A_asy(i) = A_asy_i(ti(i));
    end
    plot(ti, A_asy, 'k*');
    
    xlabel('$t$','interpreter','latex','fontsize',12);
    ylabel('Amplification','interpreter','latex','fontsize',12);
    l = legend('$A_\mathrm{num}$','$A_\mathrm{asy}$');
    set(l, 'interpreter','latex','location','best');
    
    [t_num, t_asy] = comp_times(k0, p);
    plot([t_num, t_num],[get(gca, 'ylim')], 'k--');
    plot([t_asy, t_asy],[get(gca, 'ylim')], 'k-.');

   
    figure
    plot(linspace(0,1, p.N),  ef, 'k');
    xlabel('z');
    ylabel('\check{v}_0');
    
end